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Abstract 

Suspensions with fiber-like particles in the low Reynolds number regime are modeled by two different approaches that both use 
a Lagrangian representation of individual particles. The first method is the well-established formulation based on Stokes flow 
that is formulated as integral equations. It uses a slender body approximation for the fibers to represent the interaction between 
them directly without explicitly computing the flow field. The second is a new technique using the 3D lattice Boltzmann method 
on parallel supercomputers. Here the flow computation is coupled to a computational model of the dynamics of rigid bodies 
using fluid-structure interaction techniques. Both methods can be applied to simulate fibers in fluid flow. They are carefully 
validated and compared against each other, exposing systematically their strengths and weaknesses regarding their accuracy, the 
computational cost, and possible model extensions. 
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1. Introduction 

Flows with suspended solid phase occur in many appli¬ 
cations, and thus simulation techniques for such systems 
are receiving rapidly increasing interest. An important 
and mathemactically interesting special case are fiber sus¬ 
pensions. In this article we will investigate two differ¬ 
ent computational models that represent the particulate 
solid phase in Lagrangian form. Particles will be treated 
as rigid, elongated three-dimensional geometric objects. 
Elongated is here understood as the situation that the 
shape of the particles has one dimension significantly larger 
than the others, and thus our two methods apply to sus¬ 
pensions with fibers or rods. For the fluid phase, we as¬ 
sume in this article creeping flows, i. e. the Reynolds num¬ 
ber is small and the Stokes equation can provide a suffi¬ 
ciently accurate approximation to the flow field. 
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The first method is based on a boundary integral for¬ 
mulation for Stokes flow, and approximations exploiting 
the slenderness of the suspended particles. This approach 
leads to an explicit representation of the hydrodynamic 
interactions between particles that avoids computing the 
flow field explicitly [T]. The resulting global system re¬ 
presenting the hydrodynamic interactions for this slender 
body formulation (SBF) must be solved at each time step 
during a simulation nig. 

The second method employs a full 3D Eulerian represen¬ 
tation of the fluid using the lattice Boltzmann method 
(LBM), while the particles are represented as rigid, fully 
resolved, geometric objects that can move freely through 
the simulation domain. Here fluid-structure-interaction 
(FSI) mechanisms are used to couple the flow field to the 
dynamics of the suspended particles. This approach can 
lead to very high computational cost since the mesh for 
computing the flow field must be so fine that the geome¬ 
try of the suspended objects is resolved accurately enough. 
Here using parallel supercomputers is often inevitable |3]. 

Among the many interesting effects in fiber suspensions 
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Figure 1. Tumbling motion of four spherocylinders sedimenting in Stokes flow, simulated with LBM. 


we highlight here the tumbling trajectories of sediment¬ 
ing fibers in Stokes flow, as visualized in the image se¬ 
quence in Fig. [2 Experimental results for tumbling rods 
were presented in [S]. Qualitatively correct simulations 
of this dancing motion can be achieved both with the 
LBM and the SBF models. To the best of our knowl¬ 
edge, FSI-based simulations with the LBM that show this 
phenomenon have not been reported before. 

We will find that the SBF can be more efficient and accu¬ 
rately represents the hydrodynamic interactions when the 
assumptions on which it is based are well-satisfied. The 
LBM, in contrast, can be used in more general settings 
and permits many extensions, such as simulating flows 
with higher Reynolds number and more general boundary 
conditions, and treating objects with different geometry 
than rigid fibers. The LBM based approach can finally 
also represent scenarios with particles colliding with each 
other or with the bounding walls. These possible future 
extensions justify to develop a method with higher com¬ 
putational cost. 

Besides the development of the LBM with FSI techniques 
for fiber suspensions, the primary goal of this paper is 
the systematic comparison of the two different simulation 
approaches and to assess their strengths and weaknesses. 
This comparison will also be used for the cross-validation 
of the methods. For complicated multi-physics scenarios 
such as particles in suspension, the validation increasingly 
becomes a challenge in itself, especially when analytical 
model solutions are not known, and when only sparse data 
from physical experiments exist. In this case, the compar¬ 
ison of two different simulation methodologies can serve as 
a powerful alternative to assert the correctness of the mod¬ 
els, the algorithms, and their implementation in software. 
In this paper, we study in particular the accuracy of the 
simulations by comparing the translational settling veloc¬ 
ity and the angular velocity to analytical solutions for the 
slender body motion. As a more complex scenario with 
interacting particles, we will investigate tumbling fibers. 


The simulation with four tumbling fibers in Fig.j^ illus¬ 
trates the potential of the methodology for future applica¬ 
tions that may involve many interacting particles Q This 
article, however, is restricted to the analysis and validation 
of the simpler case of two tumbling elongated particles. 

The SBF [SJ [21 [3] is an asymptotic method derived from 
the integral representation of Stokes equations and leads to 
a model consisting of a coupled system of one-dimensional 
integral equations. The model uses a discrete representa¬ 
tion of each ellipsoidal fiber in terms of the fiber center 
line and takes into account the hydrodynamic interactions 
of the fluid and the fibers. Due to the long-range nature of 
the hydrodynamic interaction, a dense system of equations 
must be solved. This method is further described in Sec.jSl 
The SBF has successfully been used for numerical simula¬ 
tions of fiber suspension in gravity induced sedimentation, 
see e.g. [HIS]. 

The LBM follows an alternative modeling paradigm that 
explicitly represents the flow field in an Eulerian way, as 
described by the Stokes or Navier-Stokes equations. A 
numerical simulation then requires the discretization on 
a grid, and the definition of interactions between parti¬ 
cles and fluid and vice versa. Such approaches have been 
studied e. g. in li HU]. Suspended particles must then 
be mapped to the grid similarly to an immersed interface 
technique m] or fictitious domain method [2]. 

In this article the LBM is employed as an alternative to a 
Navier-Stokes solver. The LBM uses a coupling between 
Eulerian fluid and Lagrangian particles via the momen¬ 
tum exchange method niinKii] and by imposing mov¬ 
ing boundary conditions on the fluid. This approach allows 
the representation of arbitrary geometric shapes, provided 
the grid resolution of the LBM method permits a suffi¬ 
ciently accurate resolution. For simulations of ensembles 
of several particles with good resolution, this leads to very 
large grids with small mesh size. These in turn lead to 


^Animation available via permalink; 
https: //wwwlO. cs .fau.de/pernialink/eehaugh3oo 
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short time steps in the LBM algorithm. Combined, these 
effects can result in very large computational cost that 
can only be provided by parallel high performance com¬ 
puting. This article therefore essentially depends on using 
the massively parallel and efficient software frameworks 
WALBerla |15j and pe [161 117] . waLBerla supports 
fluid-structure interaction with the LBM in a massively 
parallel setting HttH] by coupling it to the physics engine 
pe for rigid body dynamics. The elongated particles simu¬ 
lated in this paper are modeled as spherocylinders of fixed 
radius and with different aspect ratios. 

A validation of spherocylinder simulations with the pe in 
absence of hydrodynamic interactions was presented in 
Fischermeier et al. m- The fluid-particle interaction al¬ 
gorithm with the coupled frameworks has been developed 
and validated for charged spherical particles in microfiuid 
flows in pn] . 

The LBM was previously applied to study the sedimenta¬ 
tion of elongated rigid particles at low to moderate Reynolds 
numbers. Xia et al. m simulated the settling of sin¬ 
gle elliptical particles in a narrow channel with a two- 
dimensional multi-block LBM method to examine the wall 
influence on flow patterns for different density, aspect, 
and blocking ratios. The settling of a single spherocylin¬ 
der in a channel was studied with a two-dimensional lat¬ 
tice Boltzmann direct-forcing fictitious domain method in 
Nie et al. [22] for different solid-fluid density ratios. In [22] . 
the settling of cylindrical fibers was investigated for dif¬ 
ferent aspect ratios at moderate Reynolds numbers. The 
simulations were performed in three dimensions with the 
LB momentum-exchange method in a box moving with 
the particles. 

The LB momentum-exchange method was previously ap¬ 
plied to simulate the rotational motion of single elongated 
particles in shear flow. Ku and Lin m simulated the ro¬ 
tation of single rigid cylinder-shaped particles in planar 
Couette flow for moderate Reynolds numbers. The influ¬ 
ence of the Reynolds number and the flow confinement 
was examined in two dimensions for a given aspect ratio. 
Mao and Alexeev presented three-dimensional simula¬ 
tions of the motion of single spheroidal particles in an un¬ 
bounded shear flow at low to moderate Reynolds numbers. 
In their work, the influence of fluid and particle inertia on 
the rotational motion was investigated for different aspect 
ratios and initial orientations. 

The physical models for elongated particles in creeping 
flow are summarized in Sec.|2] Details about the SBF and 
LBM are introduced in Sec.|3]and Sec.|4| respectively. The 
simulation of a single elongated particle is validated in 
Sec.jSj and the tumbling motion of two sedimenting parti¬ 
cles is studied in Sec.jGj Finally, the findings are summa¬ 
rized in Sec. I t] and possible model extensions are outlined. 


2. Creeping flow with rod-like particles 


2.1. Fluid dynamic equations 


In this article we are concerned with the flow of small rigid 
particles suspended in a viscous incompressible fluid. The 
flow of the fluid can be described by the Navier-Stokes 
equations. 


Pf 


dVL 


= -Vp-k/i/V u-l-ft,, 

V-u = 0. 


( 1 ) 

( 2 ) 


Here, u and p denote the velocity and the pressure of the 
fluid, Pf and pf denote the dynamic viscosity and density 
of the fluid, and f;, denotes the external body force density. 

A density difference between the immersed particles and 
the surrounding fluid will give rise to a motion of the par¬ 
ticles. For particles in incompressible fluids, the gravita¬ 
tional force acting on the fluid does not need to be consid¬ 
ered explicitly. Instead, the effect of the resulting pressure 
gradient and of the gravitational force acting on the par¬ 
ticle is accounted for in the force applied to the particle 


F = (Pp - P/)gF, (3) 

where Pp is the density of the particle, V denotes its vol¬ 
ume, and g is the gravitational acceleration. 

The influence of the particle on the fluid will arise from 
the boundary condition the particle imposes on the fluid 
motion. For a rigid-body motion a no-slip condition is 
imposed meaning that the fluid velocity at the boundary of 
the particle is the same as the particle velocity. Boundary 
conditions must also be applied at the outer boundaries of 
the domain of interest. Which boundary conditions to use, 
depends on the flow case. In cases when only the motion 
of the particle is of interest, outer boundary conditions 
should be chosen that influence the flow field is as little 
as possible, e. g. periodic or free-slip boundary conditions 
(see Sec. 16.31 or Sec. l5.2|l . 


2.2. Stokes equations 


For small particles the velocity scale is usually small and 
therefore also the particle Reynolds number which is de¬ 
fined as 


Rcp = 


UL 

vf 


( 4 ) 


Here the particle Reynolds number is based on a typical 
length, L, related to the size of the particle, as well as 
on a typical velocity U of the particle, and the kinematic 
viscosity of the fluid Uf = pf /pf. In this paper, we use the 
particle diameter 2r (see Fig.[^ as a typical length scale 
for the particle Reynolds number Rcp^d- 
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When Re 1, the inertial and acceleration terms in the 
momentum equation 0 can be neglected, resulting in the 
Stokes equations 


Vp - fifV^u = fb, 

V-u = 0. 


(5) 

( 6 ) 


£.3. Single body motion at low Reynolds number 

In some simple cases, e.g. one sedimenting sphere, cylin¬ 
der or spheroid in an unbounded Stokes flow, analytical 
solutions can be found. These solutions are described be¬ 
low. More complex problems like interacting particles or 
particles in bounded fluid domains must be solved using 
computational models such as LBM or SBF. 

In the regime of low Reynolds numbers, a linear relation 
exists between the force F or the torque M applied to 
a particle and the resulting translational velocity U or 
angular velocity w, respectively. For a slender body with 
rotational symmetry around its major axis, given by a unit 
vector t, the corresponding equations are m 

U = H-iF or w = M/ 7 ^, (7) 


considering only torques applied perpendicular to t. Here, 
Jr is the rotational friction coefficient and H the transla¬ 
tional friction tensor. The latter depends on the transla¬ 
tional friction coefficients 7 ^ and for motion parallel 
and perpendicular to the symmetry axis t of the particle 
as | 2 B] 

H = 7|ltt^ + 7t^(I-tt^) , (8) 

with identity matrix I and dyadic product tt^. For forces 
applied parallel (lengthwise motion) or perpendicular (side- 
wise motion) to t, the expression for U in Eqn. 0 simpli¬ 
fies to 

Ull=F/ 7 il or U^=F/ 7 tL, (9) 


respectively. Thus, the velocity and angular velocity of 
a particle under external forces and torques can be ob¬ 
tained provided the friction coefficients are known. How¬ 
ever, these friction coefficients depend on the geometry of 
the particle and are known analytically only in a few cases 
described below. 

Cox [27] derived analytical formulas for the force acting 
on a long, slender body at rest in Stokes flow from an 
expansion of the velocity held in terms of the parameter 
E = r!L. Here, r is the radius of a particle and L its length 
(see Fig.j^. In case of a cylinder or spheroid, the transla¬ 
tional motion can be described by the friction coefficients 


7 I' = 2 ^ 


TTflfL 


ln(l/e) + Cl 


and 7 ^ = 4-" 


TTflfL 


ln(l/e) + C 2 ’ 


( 10 ) 



(a) Spherocylinder and coordinate sys¬ 
tem 


L 





(b) Cylinder 


L 





(c) Spheroid 


Figure 2. Slender body geometries and used coordinate system. 


C 2 = —l/ 2 -|-ln 2 holds, for a spheroid Ci = — 1/2 and C 2 = 
-1-1/2. The above formulas by Cox are valid for large aspect 
ratios, due to errors in the order of O ^In (e)~^^ [27| . 

In their work on diffusion of cylinders |^, Tirado et al. 
give the diffusion coefficients D for rotational diffusion, as 
well as diffusion parallel and perpendicular to the cylin¬ 
der axis. These diffusion coefficients are at given tem¬ 
perature T linked to the friction coefficients 7 via the 
Einstein-Smoluchowski relation 7 = ksT/D with ks as 
Boltzmann’s constant. Using this relation, the analytical 
derivation in leads to the translational frictional coef¬ 
ficients for lengthwise and sidewise motion, respectively 


7i' =2 


TTflfL 


In a - 


and 


7 .^ = 4 


in ( 


_L ’ 


( 11 ) 


with a = ^ = ^ • For rotational motion perpendicular to 
the cylinder axis, according to Tirado it holds that 


1 TTflfL^ 

3 In a -|- (5-*- 


( 12 ) 


The results of Tirado et al. closely resemble those of Cox, 
however, the shape-dependent corrections v and S are func¬ 
tions of a. In 1^, Tirado et al. presented the following 
relations as simple quadratic fits in 1 to the numerical 
values obtained in previous work [^ 

u'l = -0.207 -k 0.980/a - 0.133/a^ (13) 

= 0.839-k 0.185/a-k0.233/a2, (14) 

J-L = _o.662-h 0.917/a- 0.050/a2. (15) 


where the constants Ci and C 2 depend on the shape of With these corrections, the Tirado results are valid for 

the particles. For a circular cylinder Ci = -3/2-1-In 2 and aspect ratios in the range of 4 < 1/e < 60 |2H] . 
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For the friction coefficients of the spherocylinders used as 
model particles in our LBM simulation, no analytical for¬ 
mulas are known. However, the above relations for cylin¬ 
ders give a good approximation for what to expect in case 
of spherocylinders (see Sec. [5|. In this article, we will use 
the more recent results of Tirado et al. as reference. 


3. The slender body formulation 

Stokes equations can be reformulated in terms of a bound¬ 
ary integral equation. In this framework the fluid velocity 
due to the motion of a single particle can be computed by 
solving an integral equation stated solely over the particle 
surface. 

Consider a straight, rigid and slender particle of length L 
and radius r. If the particle has a large aspect ratio, i.e. 
L r, it can be referred to as a slender body. For slender 
bodies, a slender body approximation can be used. 


3.1. Non-local slender body approximation 

The slender body approximation is an asymptotic model 
derived from an integral representation of Stokes equa¬ 
tions. The model relates the velocity of the slender body’s 
surface to forces that are consistent with that motion and 
are exerted along its centerline. In the derivation, higher- 
order terms in the slenderness parameter e = rjL have 
been neglected and the accuracy of the final equation for 
the velocity of the fiber center-line is of order 0(£^lne) 
For several interacting fibers, the accuracy is of 0{e). For 
details on the derivation, see the work of Batchelor [B], 
Keller and Rubinow m, Johnson [5], and Gotz [5^ . 


3.1.1. Fiber velocities and force distribution 


Assume that we have a system of M fibers. Let the center- 
line of each fiber be parameterized by s € [—/, 1] where I is 
the half length of the fiber. For fiber m the coordinates of 
the center-line is given by Xm(t) = Xm(t) + stm(t) where 
Xm is the center point and the unit tangent vector of 
the fiber and m = 1, 2,..., M. 

Assuming that the fluid exerts a force per unit length on 
fiber m, the slender body approximation for the velocity 
of the center-line of fiber m is given by 

87r^/(xm + stm) = [d (I -t- -f 2 (I - tmt^)] fm(s) 

(16) 

+ (I + t,,OK[f,,] (s)-pV,,(s). 


Here d is a geometry parameter 


d= -ln{e'^e), 

and K [f] (s) is an integral operator given by 


K [f] (s) 




(17) 


(18) 


The contribution to the velocity of fiber m from the hy¬ 
drodynamic interaction of all other fibers in the system is 
accounted for in Vm(s) as 

M A 

= G(R|,,(s,s'))fi(s')ds'. (19) 

where Rim(s,s') = x^ + st^ — (xi -k s't|) is the distance 
between one point on fiber m and one point on fiber I. 

In free-space (no outer boundary conditions) the Green’s 
function reads 

G(R) = I Si'll+^'><'1) I;;;; ( 20 ) 


Here, S(R) = (I -k RR'^)/|R| with R = R/|R| is the 
free space Stokeslet and D(R) = (I — 3RR^)/|Rp is the 
dipole doublet. 


The unknowns in Eqn. (16) are the translational and the 
rotational velocities, x^ and tm, and the force distribution 
along the fiber fm(s). To close the formulation (161, we use 
the additional conditions stating that the integrated force 
and torque on each fiber must balance the external forces 
and torques applied to the fibers 


F = 


fm(s)ds, ^ = y s(tm X fm(s))ds. 


( 21 ) 


To solve Eqns. (16) and (21), the force on each fiber is 


expanded as a sum of Legendre polynomials. 


1 ^ 

= 2 ^ ^ ^ 

n—1 


a”Pn(s), 


( 22 ) 


where P„ is a Legendre polynomial of degree n and the co¬ 
efficients are unknown vectors with three components, 
one for each direction in space. The choice of N will be a 
parameter in the numerical method. With this approach, 
the coefficients, a”, will be given as the solution to a dense 
linear system of equations with 3MN unknowns. The sys¬ 


tem of equations is derived from Eqn. (161 using the force 
expansion Eqn. (22), orthogonality properties of Legendre 


polynomials, and the fact that the operator K in Eqn. (181 
diagonalizes under the Legendre polynomials, cf. 
details on the derivation, see [3]. 


For 


Once the system of equations for the Legendre coefficients 
has been solved, the force on each fiber can be computed, 
and the translational and rotational velocities for each 
fiber can be computed using 


' 1 /!,V4s)ds, 


S'KfJ.f 

3d 


(24) 


~^ 2 TrpfL^ ~ tmt^) sVm(s) ds. 

By integrating Eqns. (23) and ( p4| ) in time, the position 
and orientations of the fibers can be updated. 
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3.1.2. Numerical algorithm 

The numerical algorithm developed to solve this problem 
is presented in detail in |3] where also the accuracy of the 
numerical method is carefully studied. Here we will only 
give a short summary of the numerical algorithm. 

In the numerical treatment of this problem, integrals of 
the form 



G(Il{s,s'))Pkis')ds' 


Pn{s)ds 


(25) 


must be computed. Note that in order to use the Legendre 
expansion, the equations are solved in a dimensionless form 
such that — 1 < s < 1. For the inner integral in Eqn. (251 
formulas for analytic integration have been developed |3]. 
The outer integral in Eqn. (25) is evaluated numerically 
by splitting the integration interval into Nq sub intervals, 
using a three-point Gauss quadrature rule on each interval. 


The linear system of equations for the coefficients in the 
Legendre expansion is a dense system and is solved itera¬ 
tively using GMRES which on average converges (depend¬ 
ing on the distance between the fibers) within four or five 
iterations. 


To update the position of the fibers, Eqn. (231 and Eqn. (241 
are discretized in time using an explicit second-order time¬ 
stepping scheme with a fixed time step. 


3.1.3. Extension to periodic boundary conditions 


To perform simulations in a periodic domain, we must 
work with a periodized version of the Green’s function in 
Eqn. (19). This term will now also include the contribu¬ 
tion from all periodic images of the fibers. It has no closed 
analytical form, but can be thought of as a sum over an 
infinite periodic array of free space Stokeslets. 


The periodic Stokeslet is evaluated through sums in real 
and Fourier space, and although rapidly converging, it is 
more costly to evaluate than a Green’s function with a 
closed analytical expression. Therefore, parts of the peri¬ 
odic Stokeslet is initially evaluated on a uniform grid with 
grid size hg, covering the domain [Q,Lx/2] x [0, Ly/2] x 
[Q,Lz/2]. Due to symmetries, trilinear interpolation can 
be used to obtain the values needed for any coordinate 
in a periodic box of size \—L^/2^Lj;/2] x \—Ly/2,Ly/2] x 
[-L,/ 2 ,L,/ 2 ]. 

For conditions on convergence of this sum and its practical 
evaluation, see [3]. 

The integral over the periodized Green’s function is treated 
numerically in the same way as the outer integral in Eqn. (|25|). 


4. Fluid-particle interaction with the lattice Boltz¬ 
mann method 


The lattice Boltzmann method 


The lattice Boltzmann method (LBM) is a numerical scheme 
for simulating hydrodynamics based on kinetic theory of 
gases. In the LBM, the phase space is discretized into a 
Gartesian lattice fldx C of dimension D with spacing 
dx, and a finite set of Q discrete velocities Cq G R-^,g G 
{!,..., Q}. Associated with each Cq is a particle distribu¬ 
tion function (PDF) fq : X Tdt t K that represents the 

probability of an ensemble of molecules located at the lat¬ 
tice site Xi G fldx to move at that velocity. These velocities 
are chosen such that within a time increment dt = 
with discrete time Tdt = {tn : n = 0,1, 2,...} C Kg , PDFs 
can move to neighboring lattice sites, or rest at a site. 

The generalized discrete lattice Boltzmann equation [331 
134) with collision matrix S 


fqip^i “f Ogdt, fyi -)- dt) fq{pC^^ t^,) ^ (^fj 

3 


fD (26) 


is an approximation of the Boltzmann equation in dis¬ 
crete phase space that can be derived from a forward- 
difference discretization in time and a spatial upwind dis¬ 
cretization |33135] . This equation describes the advection 
of PDFs between adjacent lattice sites and subsequent col¬ 
lisions that are represented by the collision operator in the 
right hand side. The equilibrium distribution function, 
fq’^, as proposed by He and Luo m is given by 


/r(^/’") = ^9 Pf+Po 


2cf 


(c»2- 1 (u^u) 


(27) 


and recovers the incompressible Navier-Stokes (momen¬ 
tum) equation up to an error of order O(Ma^). Here, 
Ma = ^ is the Mach number, a measure for compress¬ 
ibility that depends on the characteristic velocity U of the 
fluid and the thermodynamic speed of sound Cg of the lat¬ 
tice model. The equilibrium distribution function depends 
on the local fluid density p/(xi,t„) = po + 5p(pCi,tn) with 
average value po and fluctuation 5p, and up to quadratic 
order on the local velocity u. These macroscopic quantities 
can be computed as moments of fq as 


Pf{yii,t) = Y,fq{yi„t), 


Q 

u(Xi,t) = j^J2Cqfq{Ki,t), 
q 


(28) 


and the pressure p is given by the equation of state for an 
ideal gas, p(xj,t) = c^p/(xi,t). 

We use the D3Q19 model of [3H] with c = dx/dt. 


Cs = c/v^ 


(29) 


and the following weights Wq for the different directions: 
wi = 1/3, W 2,...,7 = 1/18, and iy 8 ,...,i 9 = 1/36. 
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We employ the stable and accurate two-relaxation-time 
(TRT) collision operator by Ginzburg [5^HO] , 

E {fj - /D = {r, - /r’1+(/,“ - /r’°) • 

(30) 

with two relaxation parameters, Ae = —t ^ for even- and 
Ao for odd-order non-conserved moments. Here, Ao is a 
free parameter, and the dimensionless relaxation time t 
(or collision frequency a; = r“^) is related to the kinematic 
viscosity of the fluid by 

= (^r-c^dt. (31) 

With this definition, the LBM is second order accurate in 
space and time [33 • 

For the TRT operator, the PDFs are decomposed as 
fq = fq + fq into even and odd components 

rq = \{h + k) and /eq.« = i(/eq + /|a), 
k = hik-k) and /eq.° = i(/|q_/|a), 

with opposite velocities Cq := —Cq. The local equilibrium 
distribution function for the incompressible LBM accord¬ 
ing to m is then given for each lattice site by 



In each time step G Tdt the LBM performs a collide step 
and a stream step 

fq{l^^,t„) = fq{l^^,tn) +Xe[fqi7ii,tn) - . 

-kAo[/°(xi,t„) -/®a’°(x,,t„)] 

fqi^i T ^q^ in T dt) — _/g(Xj, tn)? (^^) 

where fq denotes the post-collision state. 

Boundary conditions are treated in the stream step by 
modifiying the post-collision states of PDFs of fluid lattice 
sites xp next to a boundary. The adaption to the bound¬ 
ary condition is performed for the PDFs associated with 
directions e^, in which the neighboring cell at x„ = x^-|-e,j 
lies on the boundary. For the simulations presented in this 
article, we apply no-slip and free-slip conditions at the do¬ 
main boundary. These boundary conditions originate from 
lattice gas bounce-back conditions and specular reflection 
conditions, respectively m- Moreover, periodic boundary 
conditions are used that cyclically extend the domain in a 
dimension. 

The no-slip boundary condition enforces zero velocity at 
the interface of two lattice cells by reverting the PDFs of 
the relevant directions as 

(36) 


For the BGK model, the effective wall locations depend on 
T. The free TRT parameter Ao allows to fix walls aligned 
with the lattice dimensions half-way between two lattice 
sites for Ao = —8(2 — w)/(8 — w) |32]. This parameter is 
used for all simulations performed in this article. The free- 
slip boundary condition enforces zero velocity in normal 
direction at the boundary, while retaining the tangential 
velocity components by reflecting the PDFs as 

/refl(g) (^F -k e, - n (n^e,) , tn -k dt) = /^(xf, t^), (37) 

where n denotes a normalized wall surface normal vector. 
The post-reflection direction associated with frefi(q) can be 
computed as ei.efi(g) =6^ — 2 (n^e^) n. 

Periodic boundary conditions are realized by streaming the 
PDFs of the relevant directions to fluid lattice sites located 
next to the boundary at the other side of the periodic 
domain. 

In a domain with periodicity in all coordinate directions, 
and when a constant force is applied to the embedded par¬ 
ticles at each time step (see Sec.[6l), the system must be 
stabilized to prevent it from accelerating infinitely. To 
keep the net momentum in the system constant, we apply 
a momentum stabilization technique: The average veloc¬ 
ity in the whole domain is computed and then subtracted 
from the macroscopic velocity when computing the equi¬ 
librium distribution function. Thus, the LBM performs 
relaxation towards a state with zero net momentum, while 
preserving all other properties. 

1-2. The momentum exchange approach 

The hydrodyamic interactions of the particles via the fluid 
with the LBM are modeled by means of the momentum 
exchange approach, as introduced by Ladd UliS]. This 
exploits the mesoscopic origin of the LBM to compute 
the momentum exchange between the fluid and the sus¬ 
pended particles directly from PDFs adjacent to the par¬ 
ticle boundary. While the original method represents par¬ 
ticles as fluid-filled shells, we use the more stable variant 
with solid particles according to Nguyen and Ladd [TS] . 

The solid, rigid objects are mapped onto the lattice such 
that the lattice sites are divided into a set b of moving 
obstacle cells whose center is overlapped by particles and a 
disjoint set F of fluid cells. Thus, particles are represented 
as obstacle sites with a staircase approximation of their 
surface. The obstacle cells b that are adjacent to a fluid 
cell F in any direction q are denoted as surface cells s. The 
momentum transfer from particles to the fluid is modeled 
by the velocity bounce-back boundary condition [T^] 

fq (xF,t„ -kdt) = fq (xF,t„) - 2 ‘^poCqUs (38) 

that adapts the fluid velocity at a given fluid cell F ad¬ 
jacent to the moving boundary to the local velocity at 


fq{xF,tn +dt) = fq{xF,tn)- 
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an obstacle surface cell s. The PDF that is bounced back 
from the obstacle is updated such that the usual no-slip 


boundary condition Eqn. (36) is recovered in the stationary 


case, i. e., when the fluid velocity matches the boundary 
velocity m 

The amount of momentum Spg transferred from the fluid 
to a particle within a time step along a given link in direc¬ 
tion Cg, can be computed (cf. [H]) as 

(5p, = [Cgfg (XF, t„) - Cgfg {x F , -f dt)] dX^ (39) 

from the difference of the momentum densities associated 
with the incoming PDF of a fluid cell F, and the PDF re¬ 
flected from the particle surface in direction Cg. The corre¬ 
sponding force acting on the particle along a given link can 
then be obtained from the relation F = ^ and Eqn. (391 


with fg given by Eqn. (381. Summing up the force con¬ 


tributions of the momenta transferred from fluid cells F 
to neighbouring surface cells s of a given particle results 
in the overall hydrodynamic force on the particle as given 
in [13] 


= EE 

s (jeD, 


2fg (xF,t„) - 2^/9 oC^Us 


(40) 


Here, xp = Xg -t-e^ and Dg is the set of direction indices q, 
in which a given s is accessed from adjacent F. The overall 


torque can be computed analogously to Eqn. (401 by 


replacing Cg by Cg x (xs — X( 7 ), with the particle’s center 
of mass xc- 


For simplicity, we use the mean density po in Eqn. (381 


and consequently Eqn. (401 instead of the fluid density pj 
in the neighbouring fluid cell. This is a good approxima¬ 
tion for incompressible LBM in absence of large pressure 
gradients. Moreover, as analysed in m, even large devia¬ 
tions from pq would have negligible effect on the accuracy 
of the hydrodynamic force computation. 


The solid particles lead to fluid cells appearing and disap¬ 
pearing due to particle movement. Lattice sites that are 
uncovered by a particle that is moving away are re-filled 
by setting the PDF at this site to the equilibrium distri¬ 
bution according to Eqn. ( [M| . We use the mean density, 
together with the particle surface velocity at that cell from 
the previous time step to compute {po, Us(xs(t„ — dt)). 


4-3. Coupling the LBM to rigid body dynamics 

For parallel fluid-particle interaction simulations, the pre¬ 
viously described methods are implemented in the LBM- 
based flow solver waLBerla that is coupled to the physics 
engine pe. Both software frameworks are designed for mas¬ 
sively parallel simulations, using a domain partitioning ap¬ 
proach for distributed memory parallelization with MPI. 
The coupling strategy and the implementation of the fluid- 
particle interaction algorithm are described in |35J US] and 
are only outlined below. 


WaLBerla |T71 [T51I20] is a parallel software framework 
for simulating fluid flow that employs the LBM. For the 
distributed memory parallelization, the simulation domain 
is decomposed into a cartesian grid of equally sized blocks 
that are assigned to the different MPI processes. On each 
process, data from adjacent lattice sites on a neighboring 
process is accessible via ghost layers. WaLBerla per¬ 
forms in each time step the streaming (Eqn. ( |3^ ) and colli¬ 
sion (Eqn. (HI) of the LBM that are fused to a performance- 
optimized stream-collide step. Incorporated in this step is 
the treatment of the boundary conditions applied at the 
domain boundary, together with the velocity bounce-back 
conditions (Eqn. (38)) at the particle surface that model 
the momentum transfer to the fluid based on the local 
particle velocities 


The pe [IHl HZ] is a framework for large-scale parallel rigid 
multi-body dynamics simulations. The rigid objects are 
geometrically fully resolved, and their translational and ro¬ 
tational motion is computed including frictional collisions 
of individual particles. Of the algorithms for multi-contact 
problems available in the pe, we employ the parallel fast 
frictional dynamics (FFD) algorithm [T6] that is based on 
Kaufman et al. jlS]. For the parallelization, the domain 
is partitioned exactly as for WaLBerla. Each pe process 
handles the particles whose centers of mass are located in 
the associated subdomain. For the collision handling, each 
process additionally stores shadow copies of intersecting 
particles [33]. A detailed description of the parallel FFD’s 
time-stepping procedure including MPI communication is 
provided in HSUS]- 


The momentum transfer to the particles is modeled by 
computing the contributions to the hydrodynamic force 
at the particle surface in WaLBerla as Eqn. (40), from 
which the pe aggregates the total force acting on the cen¬ 
ter of mass and the corresponding torque. The new posi¬ 
tions and orientations of the particles are computed in the 
subsequent pe step, together with their translational and 
angular velocities. These velocities affect the fluid motion 
in the next time step, which in turn influences the parti¬ 
cles. This interaction modelling corresponds to a two-way 
coupling. 


4.4- Parallel high performance computing for the LBM 

The explicit time discretization of the LBM restricts the 
time increment, and thus often many time steps are re¬ 
quired to simulate physically relevant phenomena. Due 
to its strictly local memory access pattern in the steam- 
collide step that involves only adjacent sites, the LBM 
allows for highly parallel simulations with excellent scala¬ 
bility. The parallel scalability of the fluid-particle interac¬ 
tion algorithm implemented in waLBerla was presented 
in |3j on up to 294 912 parallel processes. 

The LBM simulations for this article were performed on 












9 


D. Bartuschat et al. / Computers & Fluids 00 (2015) 


the high performance clusters LiMa|^of the computing cen¬ 
ter RRZE in Erlangen (Germany) and SuperMUCj^of the 
Leibniz Supercomputing Centre LRZ in Garching (Ger¬ 
many). LiMa comprises 500 compute nodes, each con¬ 
taining two Xeon 5650 ‘Westmere’ hexa-core processors 
running at 2.66 GHz and with 24 GB DDRS RAM. Super- 
MUC comprises 18 thin islands with 512 compute nodes, 
each node containing two Xeon E5-2680 ‘Sandy Bridge- 
EP’ octa-core processors that are running at 2.5 GHz and 
that have 32 GB DDRS RAM. The nodes of both parallel 
clusters are connected by a high-speed InfiniBand inter¬ 
connect . 

The technical data of the parallel LBM simulations on 
LiMa and SuperMUC for the single particle motion val¬ 
idation in Sec. [ 5 ] are summarized in Tab.[TJ These simula¬ 
tions include validations of the translational and rotational 
motion, examinations of wall effects, and flow field visual¬ 
izations. As overview, the minimal and maximal problem 
sizes and the associated parallel processes are shown, to¬ 
gether with the minimal and maximal time step numbers 
and the runtimes. 


Table 1. Parallel run data for single particle motion LBM sim¬ 
ulations. Lists cluster, problem size {Lx x Ly x Lz), numbers of 
processes (^^proc.) and time steps (t^TS), and runtime (RT) for 
validations of translational (a) and rotational (b) motion, examining 


wall influence on translational (c) and rotational (d) motion 
flow field visualization (e). 

and 


Gluster 

Lx X Ly X L 2 , 

#proc. 

#TS 

RT 



[dx] 



[h] 


Super- 

2560^ X 2688 

16 X 16 X 32 

37 000 

8.0 

& 

MUG 



88 800 

19 


LiMa 

816^ 

8 X 8 X 12 

40 000 

2.6 

D 




140000 

9.4 


LiMa 

160^ X 1200 

4 X 4 X 12 

61540 

0.9 



640^ X 1200 

8 X 8 X 12 


4.0 

C 

Super- 

832^ X 1200 



6.1 


MUG 

1920^ X 2240 

16 X 16 X 32 

86 156 

8.5 



2560^ X 2688 


70 400 

15 


LiMa 

168^ 

2x3x4 

80 000 

1.3 

a 


12963 

8 X 8 X 12 

84 000 

19 

e 

LiMa 

832^ X 1200 

8 X 8 X 12 

70 400 

7.7 


The LBM simulations of the tumbling particles in Sec.|6]for 
the domain size of 576^ lattice sites and 600 000 time steps 
were performed on LiMa on 8 x 12 x 8 processes and took 
about 16 h. The runs for the domain size of 768^ sites with 
605 000 time steps were performed on SuperMUC within 
48 h. These simulations with periodic boundary conditions 
are slower than the single particle simulations, due to the 


^WWW.rrze.fau.de/dienste/arbeiten-rechnen/hpc/systeme/ 
’'WWW. lrz.de/services/compute/supermuc/ 


average velocity computation for momentum stabilization 
that requires global MPI communication. 


5. Validation and comparison of the different meth¬ 
ods for single particle motion 

This section presents a systematic validation of the models, 
algorithms and the software used. To this end, we describe 
our findings regarding the motion of a single particle un¬ 
der constant force and torque, respectively, and compare 
the results obtained from LBM simulations to analytical 
models for slender bodies in a fluid. 


5.1. Analytical formulas 


From the analytical formulas for the motion of a slender 
body in free-space given in Sec. 12.31 theoretical values for 
the translational and angular velocities of cylinders and 
ellipsoids are computed. We compare these values to our 
simulation results for validation of the implementation and 
used methods. 


For a single fiber in a free-space setting the SBF yields 
explicit formulas for the translational and angular velocity 
of the fiber. Under a constant force F or torque M, Eqns. 
(23) and (241 simplify to 


U = x=-T^[d(I + tt^)+2(I-tt^)]F, (41) 

OTTflfL 

3^ 

a; = txt = tx -- p 7 (M X t), (42) 

ZTTfJ.fL'^ 


noting that V^, = 0 when there is only one fiber present 
in the system. These expressions are valid for an ellip¬ 
soidal particle of length L and aspect ratio 1 /e represented 
by the geometry parameter d (see Eqn. 0) moving in a 
fluid with dynamic viscosity fi /. Gonsidering a force either 
parallel or perpendicular to t, Eqn. (|4T|) gives 


U 


II _ 

SBF ~ 


2d ^ 

-F 

SiTfifL 


and 


tt-L — 

^SBF — 


rf +2 F 

8TTfJ,fL 


(43) 


These are the same expressions as derived by Gox for the 
lengthwise and sidewise motion of a spheroid, see Eqn. (101. 
For M perpendicular to t Eqn. (421 simplifies to 


WSBF = 


3d 

2tt 


M . 


(44) 


The friction coefficients for the spherocylinders that we 
use in the lattice Boltzmann simulation are not known 
analytically. However, as a first approximation we can 
compare our data to the analytical results for cylinders 
and spheroids presented in Sec. 12.31 

The results for the spherocylinders are expected to lie be¬ 
tween the results for cylinders of the same radius, but with 
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the length of the spherocylinder without the spherical end- 
caps LnC {LnC = L — 2r, see Fig. |2(a)[ ) and with the full 
spherocylinder length L. 

From the analytical expressions for cylinders derived by 
Tirado et ah, the terminal translational velocity for length¬ 
wise motion can be computed from the friction coefficients 


in Eqn. (Ill with shape-dependent correction factor in 


Eqn. (13). In case of cylinders of spherocylinder length 


without the spherical end-caps, it is denoted by ClxirnC- 
For cylinders of the full spherocylinder length, the trans¬ 
lational velocity is denoted by 

Analogously, the terminal translational velocity for side¬ 


wise motion of cylinders computed from Eqn. (Ill with 


shape-dependent correction factor in Eqn. (141 is denoted 
by for the spherocylinder length without end-caps, 

and t/rjirjj, for the full spherocylinder length. The angu¬ 
lar velocities computed from Eqn. ([^ and Eqn. (121 with 
the shape-dependent correction factor in Eqn. (15) is de¬ 


noted by WTir.nC and wxir.wC for the spherocylinder length 
without end-caps and for the full spherocylinder length, 
respectively. 


5.2. Models and parameters 

The parameters used for the validation experiments are 
chosen based on LBM requirements in terms of spatial 
and temporal resolution. To sufficiently resolve the par¬ 
ticles, the radius of the spherocylinders is kept constant 
at r = 4dx, with spatial discretization dx = 10 • 10“®m. 
The aspect ratios 1 /e are varied from 4 to 14, correspond¬ 
ing to particle lengths L of 16 dx to 56 dx. As a fluid, we 
choose water at room temperature with kinematic viscos¬ 
ity Vf = I ■ 10“®m^/s and density pf = 1 ■ lO^kg/m^, cor¬ 
responding to the dynamic viscosity Pf = 10“^kg/(ms). 
For the LBM r = 6 is chosen, which results in the time 
increment dt = 183 • 10“®s by its relation to Vf and dx 
given in Eqn. (311. 


The spherocylinders are modeled with the density Pp = 
1195 kg/m^. To all these spherocylinders with different 
aspect ratios the same force = 5.128 • 10~^°kgm/s^ is 
applied for the translational velocity validation. This force 
corresponds to the gravitational force acting on a sphere 
with r = 4dx for the given density difference between 
spherocylinders and fluid (see Eqn. (|^). The parameters 
are chosen such that all simulations are performed with 
small enough Reynolds number (Rep^d < 0.040) so that we 
are in the Stokes regime. The obtained particle Reynolds 
numbers defined in Eqn. for the terminal sedimenta¬ 
tion velocity of the spherocylinders are in the range of 
Rep,d = 0.016 ... 0.040 (see Tab. lA.7|Appendix A ). Higher 
aspect ratios than 1/e = 14 are difficult to simulate with 
the LBM at very low Reynolds numbers, due to large do¬ 
main sizes required to reduce wall effects (see Sec. 15.411 and 
decreasing dt required with increasing particle lengths (see 
Eqns. (|^, (311 and (29)). 


For the rotational motion validation the torques in Tab.|2] 
acting in x-direction are applied, leading to a rotation of 
the spherocylinder axis in the x-z plane. These torques 
are chosen such that the theoretical tip velocity Utip = 
1 • 10“® m/s is the same for all aspect ratios. This velocity 
is computed as Unp = a;Tir.wcA/2, i.e. from the angu¬ 
lar velocity according to Eqn. (121 for a cylinder with the 
spherocylinder length including the end-caps. The corre¬ 
sponding theoretical particle Reynolds number is Rep^d = 
0.009. The Reynolds numbers obtained in the simulations 
in this case are Rep^d = 0.013... 0.010 for Xje = 4... 14 
(see Tab. lA.7ll . 


Table 2. Torques in 10 [kgm^/s^] applied to spherocylinders 

of different aspect ratios 1/e for rotational velocity validation. 


1 /e 

4 

6 

8 

10 

12 

14 

M, 

12.26 

17.86 

24.62 

32.38 

41.09 

50.68 


For the LBM simulations, the particles are placed in large 
domains to minimize the wall influence and to make the 
results comparable to free space. The translational ve¬ 
locity validation with the LBM is performed in a cuboid 
domain with free-slip boundary conditions applied at the 
walls. The size of the domain in the direction of motion is 
elongated and denoted by L^. In the other dimensions it 
has the same size, i. e., = Fy. Initially, the spherocylin¬ 

der is placed at the center w.r.t. the x- and y-dimension, 
and at a distance of z = 100 r from the top boundary of 
the domain located at z = 0. Here, the coordinate sys¬ 
tem depicted in Fig. |2(a)| is used, and a slice through the 
simulation domain center along the x-z plane is shown in 
Fig-i Under the influence of F/, the spherocylinder is 
then moving downwards along the domain centerline in 
positive z-direction. 

For validation of the rotational velocity, the spherocylinder 
is placed at the center of a cubic domain, and the torques 
Mj; from Tab.[2]are applied constantly. Since these torques 
act in tangential direction to the domain boundaries, no¬ 
slip boundary conditions are applied at the walls to prevent 
the fluid from accelerating infinitely. The same setups are 
used for studying the influence of the wall on the terminal 
motion and for the flowfield visualization. 

5.3. Single particle motion results 

We validate the LBM by comparing the spherocylinder ve¬ 
locities to analytical solutions according to the theories of 
Cox and Tirado et al. for cylinders, as well as the SBF for 
ellipsoids. Moreover, the influence of the particle shape, 
fluid inertia, and wall effects on the particle velocities is 
investigated for the different methods and theories. 

For the translational velocity validations, the simulations 
are performed for the domain size [2560 dx]^ x 2688 dx for 
37 000 to 64 400 time steps for the lengthwise moving par¬ 
ticles, and for 56 200 to 88 800 time steps for the sidewise 
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motion. A high number of time steps is required until the 
particles reach steady-state—the longer the particles are, 
the more time steps are needed. For the validation of the 
angular velocity, the LBM simulations are performed for 
the domain size [816dx]^ and 40 000 to 140 000 time steps. 

The sedimentation velocities for lengthwise and sidewise 
motion are shown in Fig.[^ and Fig.[^ respectively, for 
spherocylinders, cylinders, and ellipsoids with different as¬ 
pect ratios. The angular velocities are presented in Fig.[^ 
For the considered aspect ratios, the analytical expressions 
by Tirado et al. are the most accurate, wheras the asymp¬ 
totic expressions by Cox and the SBF are valid for very 
high aspect ratios (see Sec. 12.311 . Thus, the obtained veloc¬ 
ities are normalized by the solutions C/Tir,wC and wxir.wC 
of Tirado et al. for cylinders of the full particle length, in 
order to highlight the differences in the velocities. 



Figure 3. Normalized sedimentation velocities f//C/Tir,wC 

for lengthwise orientation w.r.t. direction of force Fz = 
5.128 • 10“^^ kgm/s^ for aspect ratios 1/e and radius r = 4dx. 
Comparison of LBM velocities [2560 dx]^ x 2688dx sized 

domain with free-slip boundaries, to free-space solutions by Tirado 
for cylinders of lengths including and excluding 

spherocylinder end-caps, and to Cox and SBF (Uggp) for 

cylinders and ellipsoids of full spherocylinder lengths, respectively. 


In the figures, terminal sedimentation velocities from the 
previously described LBM simulations for lengthwise (Cbm)’ 
sidewise (t/L^jy;), and rotational (wlbm) motion are plot¬ 
ted, normalized by the corresponding Tirado velocities. 
For comparison, the normalized velocities according to 
Tirado et al. for cylinders of the same total length (denoted 
by subscript ‘Tir,wC’), and for the spherocylinder length 
without the spherical end-caps are shown (‘Tir,nC’), to¬ 
gether with the Cox results (‘Cox, cyl’) for cylinders and 
the SBF results (‘SBF’). 

The LBM results in Fig.|^to Fig .|^are mean values of par¬ 
ticle velocities after a sufficient number of time steps so 
that steady state is reached with sufficient accuracy. Due 
to obstacle mapping effects, the particle velocities fluctu- 



1/e 


Figure 4. Normalized sedimentation velocities CZ/C/xir wC 1°*' orien¬ 
tation perpendicular to Fz^ r = 4dx, and aspect ratios 1/e. Com¬ 
parison of LBM velocities (L^lbm) [2560 dx]^ x 2688 dx domain 
with free-slip boundaries, to free-space solutions by Tirado 
Cox and SBF (U^^p). 


ate as the effective particle volume varies with the number 
of overlapped particle cells (see Sec. l4.2|l . The mean values 
^LBM ‘^LBM computed from the velocities sampled 
every 200 time steps. For the validation of translational 
velocity, the last 15% of these velocity values are consid¬ 
ered and the last 50% for the validation of the angular 
velocity that converges to steady state more quickly. 


The fluctuations are computed from the minimum and 
maximum values of the considered velocities as 5u = 
(t^max ~ ^min)/t^LBM analogously. The fluctua¬ 

tions for the translational motion are lower than for the 
rotational motion. For sidewise moving spherocylinders, 
fluctuations of (5u = 0.8% occur for all aspect ratios. For 
the lengthwise moving spherocylinders, the fluctuations 
correspond to a value of 5u = 0.8% for 1/e = 4 and de¬ 
crease with increasing particle length. For the rotational 
motion, fluctuations of S^i = 8% arise for 1/e = 4 and 
decrease to = 2.8% for the longest examined parti¬ 
cle. Exact figures of the fluctuations are given in Tab. lA.7l 
(Appendix A), together with figures of the LBM results 
and associated Reynolds numbers, of the analytical solu¬ 
tions according to Tirado et al. and of the relative devia¬ 
tions Aj-U = (17lbm ~ UTir)/UTii and (defined anal¬ 
ogously) of the LBM results from the Tirado velocities. 


The terminal velocities for translational motion of sphero¬ 
cylinders in Fig.and Fig.are close to the cylinder ve¬ 
locities according to Tirado for full spherocylinder lengths. 
For small aspect ratios 1/e < 8 the spherocylinder veloc¬ 
ities are slightly higher than f/xir.wC and are thus closer 
to those of cylinders with the length of of the spherocylin¬ 
der without the end-caps. The LBM velocity for 1/e = 4 
is by Aril = 4.6% higher than f/Tir.wC for lengthwise 
motion and by 4.1% for sidewise motion (see Tab. roil . 
For large aspect ratios, the LBM sedimentation velocities 
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are slightly lower than the analytical solutions tlxir.wC for 
cylinders. The LBM velocity for the highest aspect ratio is 
by Art/ = 1.6% lower than t/xir.wC for lengthwise motion 
and by 2.3% for sidewise motion. 

The lower spherocylinder velocities for 1/e > 10 result 
from wall effects that slow down the particles compared 
to a free-space setting. These effects are examined more 
closely in Sec. 15.41 The wall influence is higher for side- 
wise orientation, resulting in LBM results and the Tirado 
results C/xir.wC to coincide for the lower aspect ratio of 
1 /e = 8 , compared to 1 /e = 10 for lengthwise orientation. 


lated with the LBM. Here, the SBF velocities of ellipsoids 
are first higher than the velocities according to Tirado for 
1 /e > 8 and converge only very slowly towards the results 
for cylinders. For sidewise motion the differences between 
the theories are smaller, and the associated velocities con¬ 
verge faster than for the lengthwise motion, due to a lower 
dependence on the particle shape. 

For the rotational motion of elongated particles shown in 
Fig .[^the differences between the theories and methods are 
larger than for translational motion. Also the Tirado re- 


Inertial effects are included in the LBM simulations and 
result in lower sedimentation velocities compared to Stokes 
flow since the fluid resistance increases with Reynolds num¬ 
ber. For both, lengthwise and sidewise motion, the LBM 
velocities tf/gM are higher than C/xir.wC only for low as¬ 
pect ratios, i.e., for minimal Reynolds numbers based on 
the particle length. For long particles J/^bm lower than 
forir.wC) and thus inertia might play an additional role in 
the translational motion simulations. 

For both orientations, the velocities according to Cox and 
the SBF exhibit a similar behavior for small aspect ra¬ 
tios. Both velocities increase at the same rate with in¬ 
creasing particle lengths for low aspect ratios and slowly 
decrease for higher aspect ratios. However, while Cox’ 
results quickly converge to values close to the Tirado ve¬ 
locities, the SBF velocities of ellipsoids stay higher. For 
8 < 1/s < 14, the SBF velocities are by 10 — 15% higher 
than the results by Tirado. This effect can be attributed to 



lA 


Figure 5. Normalized sedimentation velocities F/LxirwC for 
lengthwise orientation w.r.t. Fz for r = 4dx and high aspect ra¬ 
tios 1/e. Comparison of theories for cylinders by Tirado 

^Tir nc) Cox and for ellipsoids by the SBF (Uggp). 

the different particle shapes, which becomes clear in Fig.[^ 
that shows the velocities according to the different theories 
for lengthwise motion at higher aspect ratios than simu¬ 



Figure 6. Normalized angular velocities a;/t<Jxir,wC of particles with 
r = 4dx and aspect ratios 1/e, resulting from torques Mj, in Tab. [2] 
Comparison of LBM velocities in [816 dx]^ domain with no¬ 

slip boundaries, to free-space solutions by Tirado (a;xir,wC) <^Tir nc) 
and the SBF (tasHp). 

suits depend more strongly on the particle length. As for 
translational motion, the velocity according to the SBF 
is higher than the Tirado solution for cylinders with full 
spherocylinder length. The angular velocities of the sphe- 
rocylinders from the LBM simulations lie between the so¬ 
lutions for cylinders of a length with and without the end- 
caps for all aspect ratios. Again, the LBM velocities are 
closer to the values of a cylinder of the full spherocylin¬ 
der length. However, the relative deviation from wxir.wC is 
significantly higher than for the translational motion. For 
1/e = 4, Ar w = 45% and for 1/e = 14 still Ar^ui = 12%. 
With increasing aspect ratio, the LBM velocities more 
closely approach the Tirado velocities with the full sphe¬ 
rocylinder length than the SBF results do. The angular 
velocities by the SBF exceed the Tirado results by almost 
90% for 1/e = 4 and still by more than 50% for 1/e = 14. 
The wall influence is negligible for rotational motion, as 
shown in Sec. 15.41 With increasing particle length, the 
angular LBM velocities w£bm decrease towards WTir.wC- 
Thus, also for rotational motion the LBM results might 
be influenced by inertial effects that lead to lower veloci¬ 
ties, in addition to the higher hydrodynamic similarity of 
longer spherocylinders and cylinders of same length. 
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5.4- Wall influence in LBM simulations 

The influence of the wall effect in the LBM simulations on 
the settling velocity of spherocylinders is examined for do¬ 
mains with quadratic cross-section of = Ly := (160, 320, 
480, 640, 832) dx and constant length = 1200 dx, and 
with sizes [1280 dx]^ x 1600 dx, [1920 dx]^ x 2240 dx, and 
[2560 dx]^ X 2688 dx. The simulations investigating the wall 
influence on rotational motion are performed for cubic do¬ 
mains with edge lengths of (168, 324,480,648, 816,1296) dx. 
For the aspect ratios 1/e = 8 and 1/e = 12, 80 000 and 
84 000 time steps are performed, respectively. 

The terminal translational and rotational velocities are 
presented in Fig.on the left and right ordinate, respec¬ 
tively, dependent on = Ly plotted as on the ab¬ 

scissa. The displayed velocities are again the mean values 
of the particle velocities evaluated every 200 time steps. 
For the translational motion, the last 34% of these ve¬ 
locities are considered, and the last 15% for the largest 
domain only. For the rotational motion, the last 50% are 
considered. 


t/[10-3 m/s] w [1/s] 



Figure 7. Sedimentation and rotation velocities of spherocylinders 
with radius r = 4dx and aspect ratio 1/e = 8—and 1/e = 12 for 
rotation—from LBM simulations in water-filled domains of different 
sizes (Lx,y)- Domains for the sedimentation velocities of length- 
wise (F]^g]y[) and perpendicular oriented spherocylinders 

(w.r.t. force F^) are cuboids with free-slip boundaries, and cubes 
with no-slip boundaries for the angular velocities resulting 

from torques Mx in Tab.[^ 


The terminal translational velocities of the spherocylinders 
with 1/e = 8 in Fig.[^ increase signiflcantly with domain 
size. For the translational motion of lengthwise oriented 
spherocylinders, the retarding effect of the confined do¬ 
main on the sedimentation velocity is negligible for do¬ 
main sizes of [1280 dx]^ x 1600dx and above. For sidewise 
oriented particles, only the [2560dx]^ x 2688dx domain is 
sufficiently large to assume negligible wall effects. These 


findings confirm the statement in Sec. l5.3l that the wall ef¬ 
fect is measurable for 1/e > 10 in the [2560 dx]^ x 2688 dx 
domain, and that the corresponding velocities are under¬ 
estimated. 

For the rotational motion, there is hardly any impact of the 
confined domains on the angular velocities for both aspect 
ratios 1/e = 8 and 1/e = 12. Thus, the wall effect for the 
[816dx]3 domain in Sec. 15.31 is negligible for all considered 
aspect ratios. According to the angular velocities shown 
in Fig.|^ domains with edge length 480 dx and 648 dx are 
sufficient for 1/e = 8 and 1/e = 12, respectively. 


5.5. Flow field around sedimenting particles 


We present the flow field around a single sedimenting sphe- 
rocylinder in a closed domain with free-slip boundary con¬ 
ditions simulated with the LBM. The simulations are per¬ 
formed for spherocylinders of aspect ratio 1/e = 12 in a 
domain of size [832dx]^ x 1200dx for 74 800 time steps. 
The fluid velocity around the sedimenting spherocylinder 
and the staircase-approximated spherocylinder itself de¬ 
picted in Fig. 1^ are visualized with ParaView. The flow 
field around the spherocylinder is shown along the x-z 
plane through the domain center for a particle oriented 
lengthwise and sidewise along the direction of the applied 
force Fz in Fig. |8(a)| and Fig. |8(b)| respectively. Since the 
velocity magnitude decays quickly with increasing distance 
from the spherocylinders, white isosurface contour lines 
indicate the velocity magnitude with logarithmic contour 
intervals. The flow direction is depicted by small streaks 
of uniform length. The sixteen contour lines are plotted 
for both orientations at velocity magnitudes in the range 
of 304 • 10“® m/s to 304 • lO-^^ m/s. 


The initially quiescent fluid is dragged along by the sphero¬ 
cylinder and starts moving adjacent to the particle. Along 
the channel centerline a fluid flow develops in the direc¬ 
tion of spherocylinder motion. The walls in z-direction 
cause a counterflux and thus a vortex. Initially the vortex 
is located next to the particle and then moves to a posi¬ 
tion half-way between particle surface and confining walls. 
For the lengthwise oriented spherocylinder in Fig. 8(a) the 
contour lines have a similar shape as the particle. For the 
sidewise oriented spherocylinder in Fig. |8(b)1 the contour 
lines next to the particle are circular, and they have simi¬ 
lar shape as for the lengthwise oriented spherocylinder in 
some distance. However, the shape is more circular, i. e., 
the flow is influenced to a higher degree also orthogonal to 
the movement direction. Moreover, the sidewise moving 
spherocylinder influences the fluid at larger distances than 
the lengthwise oriented spherocylinder, as can be seen from 
the distance of equivalent contour lines to the particles. 


In accordance with the results in Sec. 15.31 the particle with 
sidewise orientation moves slower than for lengthwise ori¬ 
entation. The same position in movement direction is 
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(a) Lengthwise moving spherocylinder after 22 000 time 
steps. 


(b) Sidewise moving spherocylinder after 28 000 time 
steps. 


Figure 8. Flow around sedimenting spherocylinders along x-z plane through domain center. Flow direction is indicated by small streaks, 
and velocity magnitude by color in range of 360 ■ 10“® m/s (red) to Om/s (blue) and logarithmic contour lines. 


reached after 22 000 time steps for lengthwise orientation, 
and after 28 000 time steps for sidewise motion. 

The flow field caused by the sidewise moving spherocylin¬ 
der is nearly equal in both directions orthogonal to the 
movement direction. Only very close to the spherocylin¬ 
der, the fluid velocity in x-z plane shown in Fig. |8(b)1 is 
found to be slightly higher than in the orthogonal y-z 
plane. 

6. Tumbling particles 

In this section we study the tumbling motion of two elon¬ 
gated particles in a periodic domain. We present simu¬ 
lation results for spherocylinders with the LBM, as well 
as SBF results for ellipsoidal fibers. A comparison of the 
results from the two methods is presented for different do¬ 
main sizes and aspect ratios. For the LBM, the influence 
of the initial particle distance on the tumbling behavior is 
examined. 

6.). Setup and parameters 

The LBM simulations are performed with two spherocylin¬ 
ders of radius r = 4dx, with dx = 4.98 • 10“® m. The par¬ 
ticle length, L, is varied between 40 dx and 56 dx, cor¬ 


responding to the aspect ratios 1/e = 10 to 1/e = 14. 
Initially the particles are placed, aligned with the direc¬ 
tion of gravity, in a periodic box at the center w.r.t. y- 
direction and centered around the middle of the domain 
in x-direction with a given center-to-center distance dist 
w.r.t. the coordinate system depicted in Fig. |2(a)| The 
dimensions of the box are 576 dx or 768 dx in all three 
directions, i.e. = Ly = L^. For the smaller domain, 

600 000 time steps are performed and for the larger domain 
605 000 time steps. To keep the system from accelerating 
infinitely due to the constantly applied force, the momen¬ 
tum stabilization technique described in Sec. l4.1l is applied. 

As fluid, we model water at room temperature, with the 
kinematic viscosity and the density Pf^H 20 from 

Sec. 15.21 The density of the particles is set to Pp = 1492 kg/m^. 
The force acting on the particles in z-direction caused by 
gravity results in the values presented in Tab.jS] 

Table 3. Gravitational forces applied to the particless of different 
aspect ratios 1/s in z-direction for the tumbling simulations. 


l/e 

10 

12 

14 

F, [10-®kgm/s2] 

1.119 

1.358 

1.598 


The resulting particle Reynolds number Rep^^; for the par- 
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tides based on the mean sedimentation velocity U* for 
the LBM lies in the range 0.052 to 0.063 for 1/e = 10 and 
1/e = 14, respectively (see Tab.|5]), so that we are essen¬ 
tially in the Stokes regime. For the LBM simulations, we 
apply the relaxation time r = 6 that results in the time 
increment dt = 4.55 • 10“^ s (see Eqn. @). 

In the SBF we use the following parameters = 5, 48 
quadrature points along the fiber (see Eqn. ([2^), and dt « 
0.003 s. 

An SBF simulation with 8300 time steps (TS) was per- 
fomed within 66.40 s on a single core of an Intel Core i7 
4770 processor running at 3.4 GHz, i. e., with 125TS/s. 
The parallel runtimes for the LBM on 768 cores given in 
Sec. 14.41 correspond to 10.4TS/s on LiMa and 3.5TS/s on 
SuperMUC. 

6.2. Flow field around tumbling spherocylinders 

We present the flow field around two tumbling spherocylin¬ 
ders simulated with the LBM in a periodic domain of size 
[576 dx]^, filled with water. The setup and parameters de¬ 
scribed in Sec. 16.II are used for spherocylinders of aspect 
ratio 1/e = 12 and an initial center-to-center distance of 
16 dx, with dx = 4.98 • 10"® m. 

In Fig.j^the flow field is visualized along a slice in the x-z 
plane through the domain center. The velocity magnitude 
is represented by a color range from red to blue and by 
(sixteen) white isosurface contour lines of logarithmic in¬ 
tervals in the range of 2.60 • 10“®m/s to 96.1 • 10“®m/s. 
The flow direction is represented by streaks of uniform 
length. In addition to the image sequence in Fig.[^ an 
animation is available via a permalinl^ 

At the beginning of the tumbling period, both particles 
are aligned with the external force direction (see Fig. |9(a)| , 
and then start rotating (see Fig. |9(b)[ ). As the particles ro¬ 
tate, they move apart in x-direction until they are oriented 
perpendicular to the force (see Fig. |9(c)[ ) and reach the 
maximum distance. The rotation continues, and the par¬ 
ticles move further along the force direction (see Fig. |9(d)] ) 
until the next period begins. A more detailed description 
of the particle positions and velocities is given in Sec. 16.31 
The fluid is dragged along with the spherocylinders, while 
the flow velocity quickly decreases with increasing distance 
from the particles. As expected, the flow field is at all 
times symmetric w.r.t. the domain center in x- and y- 
direction and changes in z-direction. During the tumbling 
motion, two vortices are forming at the same height (i. e., 
z-coordinate) as the particle center. Half-way between pe¬ 
riodically following particles in z-direction, a region of zero 
flow velocity appears at the interface of two vortices associ¬ 
ated with these periodically neighboring particles. As the 
particles change their mutual distance, the vortices and 
the zero flow velocity regions only slightly change their 
positions in x-direction. 


^https://wwwlO.cs.fau.de/permalink/eethegh4sh 


6.3. Tumbling results for spherocylinders using LBM 

We examine the dependence of the tumbling motion on 
the spherocylinder aspect ratio, the domain size, and the 
initial distance between the particles. The influence of the 
aspect ratio on the tumbling behavior is analyzed for sphe¬ 
rocylinders with 1/e = 10, 1/e = 12, and 1/e = 14. For 
one particle, the position in x-direction and its velocity in 
X- and z-direction are presented in Fig. |10| over time. The 
simulations are performed for the domain size [576 dx]® 
and the initial center-to-center particle distance of 16 dx 
in x-direction. 

Initially, the particles move away from the domain cen¬ 
ter and from each other in x-direction. Once the maximal 
distance is reached, the particles move back towards each 
other. This motion is repeated periodically. While the 
particles aligned with the direction of gravity move apart, 
the velocity in x-direction increases, and the particles start 
rotating. Before the particles are oriented perpendicular 
to the gravitational direction, the maximal velocity in x- 
direction is reached, and the particle velocity in x-direction 
decreases again. The velocity in x-direction becomes zero 
once the particles are oriented perpendicular to gravity. 
While the particles continue rotating, the velocity in x- 
direction becomes negative, i. e., the particles approach 
each other. The velocity in negative x-direction increases 
at first and then decreases again as the particles rotate 
further. Once the particles are aligned with the direction 
of gravity, the velocity in x-direction is zero again, and the 
minimal separation is reached. In z-direction, the maxi¬ 
mal velocity occurs when the particles are oriented with 
the direction of gravity and the particle distance is mini¬ 
mal. The minimal velocity in z-direction is reached when 
the particles are oriented perpendicular to the direction 
of gravity and the particle distance is maximal. In both 
cases, the velocity in x-direction is zero. The values for 
the y-direction are not shown since the particles move in 
the x-z plane, keeping their y-coordinate. 

The maximal separation in x-direction, the maximal ve¬ 
locity magnitude in this direction, and the sedimentation 
velocity in z-direction all depend on the particle length. 
With increasing particle length, the terminal maximal sep¬ 
aration distance, the maximal velocity magnitude in x- 
direction, and the sedimentation velocity Uz increase. The 
minimal separation is nearly equal for all particle lengths. 
For 1/e = 10, the maximal separation distance per pe¬ 
riod is small at first and then increases over many periods, 
while for 1/e = 12 this distance increases only slightly, 
and the terminal separation distance is reached after five 
revolutions. For 1/e = 14, the separation distance is rel¬ 
atively high in the first period and then decreases until 
the terminal distance is reached after approximately three 
revolutions. Before examining these dependencies further, 
we analyze how the parameters depend on the initial sep¬ 
aration distance between the particles and on the domain 
size. 
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(c) Flow field after 27 000 time steps (d) Flow field after 44 000 time steps 


Figure 9. Flow field around two tumbling spherocylinders with radius r = 4dx (dx = 4.98 ■ 10“® m) and initial distance 16 dx in periodic 
domain of size [576dx]®. Flow direction is indicated by streaks, and velocity magnitude logarithmic contour lines and color in range of 
3.16 ■ 10“® m/s (red) to 20 ■ 10“® m/s (blue). 


In Fig. 11 we present LBM results for the tumbling mo¬ 
tion of particles with aspect ratio Xje = 12 for a smaller 
initial center-to-center distance of 14.8 dx. Moreover, we 
compare the particle motion for the previous domain size 
[576dx]^ and a larger domain of size [768dx]^. All the 
other parameters are kept the same as for the previous 
simulations. The particle motion conforms to the previ¬ 
ously described behavior for the different aspect ratios. 
The smaller initial distance of the particles leads to a 
smaller maximal separation at the beginning of the simula¬ 
tion. Also the tumbling period time is at first smaller than 
for the higher initial distance. The velocity in x-direction 


does not depend on the initial distance, only the sedimen¬ 
tation velocity is a bit higher at first. After a few periods, 
however, the tumbling motion for the different initial sep¬ 
arations no longer differs. 

When the domain size is increased from 576 dx to 768 dx, 
the particle motion in x-direction is not affected. I. e., the 
size 576 dx suffices for negligible mutual influence of the 
periodically neighboring particles in x-direction. Only in 
z-direction the velocity is smaller for the larger domain, 
due to the periodic boundaries. 

Among the LBM tumbling simulations presented in this 
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Figure 10. LBM results for tumbling motion of two spherocylinders 
with radius r = 4dx (dx = 4.98 ■ 10~® m) and aspect ratios 1/e = 10, 
1/e = 12, and 1/e = 14 in periodic domain of size [576 dxj® filled 
with water. Shown are the position of one particle in x-direction and 
the particle velocities in x- and z-direction over time. The spatial 
discretization is dx = 4.98 ■ 10“® m, and the initial particle distance 
is 16 dx. 


section, the change in separation distance is the highest for 
1/e = 10, and for 1/e = 12 with initial distance 14.8 dx. 
For these simulations, the measured tumbling period time 
T* and period distance D* in z-direction, and the resulting 
mean sedimentation velocity U* are shown in Tab.lH For 
both aspect ratios, the tumbling period distance increases 
as the particles move apart in x-direction. The tumbling 
period time increases at almost the same rate, leading to 


Figure 11. LBM results for tumbling motion of two spherocylinders 
with radius r = 4dx and aspect ratio l/e = 12, for initial particle 
distances disti = 14.8 dx and dist 2 = 16 dx in water-filled periodic 
domains of size [576dxj® and [768dx[®. The larger domain is indi¬ 
cated by ‘L’. 


an only slightly decreasing mean sedimentation velocity 
in both cases. The terminal sedimentation velocity for 
1 /e = 12 is higher than for 1/e = 10. 

The main characteristic parameters for the tumbling mo¬ 
tion of two spherocylinders obtained from the LBM sim¬ 
ulations are summarized in Tab.[5l The terminal mean 
sedimentation velocity U* per period rises with increas¬ 
ing aspect ratio or particle length, which agrees with the 
experimental results in [^. This behavior is contrary to 
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Table 4. Changes in sedimentation distance (D*) per period, pe¬ 
riod time (T*), and mean sedimentation velocity (C*) for differ¬ 
ent periods of two tumbling particles with radius r = 4dx (dx = 
4.98- 10“®m). LBM results for particles with 1/e = 10 and ini¬ 
tial distance 16 dx, and for 1/e = 12 with initial distance 14.8 dx in 
periodic domain of size [576 dx]^. 


1/e period 

1 

2 

3 

4 

5 

6 

7 

D* [10^® m] 

4.38 

4.78 

4.97 

5.25 

5.41 

5.54 

5.72 

10 r*[s] 

2.71 

3.52 

3.69 

3.94 

4.09 

4.23 

4.35 

U* [10-3 m/s] 

1.61 

1.36 

1.35 

1.33 

1.32 

1.31 

1.31 

D* [10-3 m] 

4.18 

4.50 

4.75 

5.02 

5.24 

5.42 

5.51 

12 r*[s] 

2.13 

2.92 

3.08 

3.34 

3.53 

3.70 

3.76 

U* [10-3 m/s] 

1.96 

1.54 

1.54 

1.50 

1.49 

1.47 

1.47 


the sedimentation results in Sec.|5] where the same force 
is applied to all particles and can thus be attributed to 
the increasing gravitational force acting on the particles. 
For the larger domain, [/* increases slightly, as observed 
previously for Uz- While the terminal mean sedimenta¬ 
tion velocity is reached for all simulations, the tumbling 
period distance D* and time T* still change for 1/e = 10 
and for 1/e = 12 with dist = 14.8 dx (see also Tab.|4|). 
For the latter, the values converge towards the results for 
dist = 16 dx. The period time decreases with increasing 
aspect ratio, whereas no clear trend is recognizable for D*. 
For the larger domain and 1/e = 12, the smaller sedimen¬ 
tation velocity results mainly from an increased distance 
D*. 

The velocities in x-direction are significantly smaller than 
the sedimentation velocities. As reported above, the max¬ 
imal velocity magnitudes in x-direction increase with par¬ 
ticle length. However, these velocities do not depend on 
the domain size and the initial separation distance. The 
maximal magnitude Ux, max when the particles separate 
is higher than the magnitude Mx, min for the approaching 


Table 5. Parameters obtained from LBM simulations of two tum¬ 
bling particles with radius r = 4dx (dx = 4.98 • 10“^ m) in periodic 
cubic domain filled with water. Sedimentation distance (D*), pe¬ 
riod time (T*), and mean velocity (U*) in z-direction, minimum and 
maximum velocities in z-direction (uz, mini max) and x-direction 
min, "i^x, max) per revolution, are shown for different aspect ra¬ 
tios (1/e), domain sizes {Lx,y,z), and initial particle distances (dist). 


1 /e 


10 

12 

12 

14 

12 

^x,y.,z 

[dx] 

576 

576 

576 

576 

768 

dist 

[dx] 

16.0 

14.8 

16.0 

16.0 

16.0 

D* 

[10-3 m] 

5.72 

5.63 

5.63 

5.92 

5.71 


[s] 

4.35 

3.82 

3.86 

3.75 

3.85 

U* 

[10-3 m/s] 

1.31 

1.47 

1.46 

1.58 

1.48 

^z, min 

[10-3 m/s] 

1.11 

1.20 

1.19 

1.25 

1.23 

^z, max 

[10-3 m/s] 

1.86 

2.09 

2.07 

2.26 

2.10 

^x, min 

[10-3 m/s] 

-152 

-179 

-178 

-199 

-177 

^x, max 

[10-3 m/s] 

164 

196 

196 

222 

196 


particles. Moreover, the maximal distance between the 
particles in x-direction converges to a value that depends 
on the aspect ratio, but not on the initial separation (see 
Fig.and Fig. [TT|. The minimal distance is hardly af¬ 
fected by these parameters. 

6 . 4 . Tumhling results for ellipsoids using SBF and com¬ 
parison to LBM results 

The slender body formulation results for two tumbling par¬ 
ticles in a periodic domain are presented in Tab.for dif¬ 
ferent domain sizes and particle distances equivalent to 
the LBM simulations. Similar to the LBM results, both 


Table 6. Simulation results using the slender body formulation 
for two tumbling fibers in a periodic domain filled with water. The 
spherocylinder radius is r = 4dx. 


1 /e 


10 

12 

12 

14 

12 

Lx,y,z 

[dx] 

576 

576 

576 

576 

768 

dist 

[dx] 

15.0 

14.9 

15.1 

15.1 

15.1 

D* 

[10-3 m] 

6.40 

6.17 

7.20 

9.01 

5.31 


[s] 

4.50 

3.91 

4.53 

5.73 

3.84 

U* 

[10-3 m/s] 

1.42 

1.57 

1.53 

1.58 

1.61 

^z, min 

[10-3 m/s] 

1.23 

1.32 

1.31 

1.35 

1.36 

^z, max 

[10-3 m/s] 

2.00 

2.36 

2.35 

2.58 

2.40 

^x, min 

[10-3 m/s] 

-180 

-224 

-224 

-260 

-224 

^x, max 

[10-3 m/s] 

180 

224 

224 

260 

224 


Uz and U* increase with the aspect ratio and with the 
domain size for constant aspect ratio. However, these ve¬ 
locities are more sensitive to the domain size than for the 
LBM. Moreover, the initial distance influences U* stronger 
than for the LBM, whereas Uz is hardly affected. In con¬ 
trast to the LBM results, the SBF period time T* does not 
systematically depend on the aspect ratio whereas the pe¬ 
riod distance D* tends to increase with 1/e. However, D* 
depends strongly on both domain size and initial separa¬ 
tion distance. With increasing aspect ratio, the velocities 
Ux increase but are not sensitive to initial distances and 
domain sizes, like for the LBM. 

A comparison of the values presented in Tab. [6] except 
for T* and U* to the corresponding values in Tab. [5] re¬ 
veals a relative difference of approximately 10 — 20% for 
1/e = 10,12. This difference can be attributed to the dif¬ 
ferent shapes of the particles in the LBM and the SBF 
respectively, see the discussion in Sec. 15.31 However, for 
1/e = 14 the difference is larger, especially for quantities 
related to the periodicity of the tumbling motion, T* and 
D*. For these quantities the relative difference is now in 
the order of 35%. In addition to the effect of different 
particle shape, the combined effect of particle and fluid 
inertia could play a role here that are both only present in 
the LBM. Since the velocity of the tumbling particles in¬ 
creases with growing particle lengths, so does the Reynolds 


18 



































19 


D. Bartuschat et al. / Computers & Fluids 00 (2015) 1-123} 


number as well as the importance of fluid inertia. 
Furthermore, all values obtained with the LBM simula¬ 
tions are consistently lower than the values obtained with 
the SBF. All of the above is in line with the findings for 
one sedimenting particle presented in Sec. 15.81 


Fig.[TT] we see that the difference caused by starting with 
two particles at two different center-to-center distances de¬ 
creases with time, and the curves will eventually coincide. 
In the SBF, the difference in the resulting curves stay the 
same throughout the whole simulation. 


Another difference between the two methods that also can 
be attributed to inertial effects only being present in the 
LBM simulations, is that the SBF results depend strongly 
on the initial center-to-center distance between the parti¬ 
cles, see Fig.jl^ For the LBM the situation is different. In 




0 5 10 15 20 25 

t [s] 


Figure 12. SBF results for tumbling motion of two fibers with 
radius r = 4dx and aspect ratios 1/e = 12 for initial fiber distances 
of disti = 14.9 dx and dist 2 = 15.1 dx. The size of the periodic 
domain is [576dx]®. 


In the upper diagram of Fig.we show a comparison of 
the sedimentation velocity for a particle with 1/e = 12 
obtained with the LBM and the SBF for initial distances 
16 dx and 14.9 dx, respectively. The results agree well on a 
qualitative level. However, for the SBF the sedimentation 
velocity is larger and the period time is slightly longer than 
for the LBM. This is consistent with the results for single 
particles presented in Sec. 15.31 where the particle velocities 
of the SBF for ellipsoids are higher than those of the LBM 
for spherocylinders, due to the different particle shapes. 

6.5. Comparison of tumbling orbits from LBM and SBF 

The tumbling motion of elongated particles can also be 
characterized by orbits in phase space, i. e., the space of 
particle positions and velocities. Since the particle motion 
is periodic in the direction perpendicular to the sedimen¬ 
tation direction, the orbits in phase space are represented 
by closed curves. We present the periodic orbits for LBM 
and SBF simulations. 

The tumbling orbits obtained from the LBM simulations 
are shown in the middle diagram of Fig.j^for spherocylin¬ 
ders in a periodic, cubic domain. As shown in Fig. |10| the 
LBM simulations converge towards periodic motion in x- 
direction after a few revolutions, i. e., the orbits are asymp¬ 
totically periodic. This effect can be attributed to the iner¬ 
tia considered by the LBM. For the aspect ratio 1/ff = 10, 
the convergence takes the longest, with increasingly larger 
cycles in x-direction. The reason for this convergence be¬ 
havior is that the initial separation distance for 1/e = 10 
differs the most from the preferred distance (see Fig.[T0|). 
The next longer particles with 1/e = 12 converge faster, 
and the maximal separation in x-direction gets larger at 
first, too. For 1/e = 14, the maximal particle separation 
at the first revolution is the highest but converges quickly. 
Both, terminal maximal separation and velocity magni¬ 
tude, increase with particle length. However, minimal po¬ 
sitions and maximal velocity magnitudes in x-direction do 
not change with the revolutions. Moreover, the orbits are 
not symmetric w.r.t. the ordinate because the maximal 
separation velocity is higher than the maximal approach 
velocity. 

The periodic orbits obtained from the SBF simulations are 
shown in the lower diagram in Fig. |13| These orbits exhibit 
in principle the same behavior as the LBM results. The 
shape of the curves is similar, and the maximal separa¬ 
tion distance increases with the particle length. However, 
the initial minimum separation distance is retained and 
thus differs for the aspect ratio 1/e = 10 from the longer 
particles. Moreover, the terminal periodic orbit is already 
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obtained after the first revolution as there is no inertia, 
and the maximal separation and approach velocities have 
the same magnitude for a given particle length. 





Figure 13. Comparison of the results from LBM and SFB for 
tumbling motion of two particles with radius r = 4 dx and different 
aspect ratios 1/e and initial distance 16 dx in a periodic domain of 
size [576 dx]®. The upper diagram shows the particle velocity in z- 
direction over time for 1/e = 12. The middle and lower diagrams 
show the tumbling orbits in phase-space (x-position vs. x-velocity) 
for LBM and SBF, respectively. 


7. Conclusions 

In this article, we have validated and compared different 
methods for elongated particles in creeping flows. The 
SBF models hydrodynamic interactions of elongated ellip¬ 
soidal particles in flow based on an asymptotic formulation 
in the slenderness parameter. Hence the model accuracy 
increases with aspect ratio. The LBM simulates hydrody¬ 
namic interactions by an explicit computation of the flow 
field and models the momentum exchange between fluid 
and particles. Here, the flow is represented on a lattice, 
and the particles are modeled as spherocylinders. 

To validate the LBM simulations, LBM and SBF results 
for the translational and rotational motion of single elon¬ 
gated particles were compared to analytical solutions for 
cylinders according to the theories of Cox m and Tirado 
et al. [2H] • While the different theories for cylinders and the 
results for spherocylinders converge quickly with increas¬ 
ing aspect ratio 1/e, the SBF solutions for ellipsoids ap¬ 
proach these results only for extremely high aspect ratios. 
The particle shape is identified as the dominant factor for 
the different velocities of (sphero-)cylinders and ellipsoids. 
This shape effect explains the consistently lower velocities 
for the LBM than for the SBF that deviate by about 10% 
to 16% for translation and by about 12% to 13% for ro¬ 
tation at 1/e > 8. Overall, the terminal velocities of the 
different theories converge faster for sidewise orientation 
w.r.t. the applied force than for lengthwise orientation. 
The translational velocities for spherocylinders agree very 
well with the theory of Tirado et al. for cylinders of full 
spherocylinder lengths, apart from retarding wall effects 
in the LBM simulations at high aspect ratios. For rota¬ 
tional motion, the angular velocities of spherocylinders are 
slightly higher than predicted by Tirado et al. for cylin¬ 
ders of full spherocylinder length. We find that in the LBM 
simulations of single particles, inertial effects play only a 
minor role compared to shape and wall effects. Finally, 
the flow field around sidewise and lengthwise moving sphe¬ 
rocylinders obtained from LBM simulations is compared. 
The comparison shows that for sidewise motion the fluid 
is affected at higher distances, leading to lower sedimen¬ 
tation velocities. Moreover, the counterflux of fluid in the 
closed domain becomes visible that results in the retarding 
wall effect. 

After the single particle motion validation, the tumbling 
motion of two elongated particles was examined in a pe¬ 
riodic domain. For the LBM simulations, the flow field 
around the tumbling particles was visualized. The parti¬ 
cle trajectories obtained from both, LBM and SBF simula¬ 
tions, were compared in diagrams, including a representa¬ 
tion of the tumbling orbits in phase-space. The qualitative 
agreement of the LBM and SBF simulation results is good, 
whereas quantitative differences become apparent. Con¬ 
sistent with the single particle motion, the characteristic 
tumbling parameters differ by about 10 — 20% for aspect 
ratios of 1/e = 10,12. For an aspect ratio of 14 the dif- 
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ferences increase, which might partially be attributed to 
higher differences in the velocities of single lengthwise sed¬ 
imenting particles. 

In contrast to the SBF, a preferred maximal particle sep¬ 
aration distance exists for the LBM in the direction per¬ 
pendicular to gravity, which depends on the aspect ratio 
and can be attributed to inertia. This effect is similar to 
the drift to stable orbits due to inertia examined in Mao & 
Alexeev m] for the motion of spheroidal particles in shear 
flow. Moreover, Jung et al. [5] report that the experimen¬ 
tal tumbling results for rods are insensitve to the initial 
separation and converge to the preferred distance. Also 
the increasing mean sedimentation velocities with particle 
length observed in the LBM and SBF simulations agree 
with the experimental investigations in [5]. 

The strengths of the SBF compared to the LBM lie in 
its significantly lower computational effort as well as in 
its ability to simulate zero Reynolds number flows. How¬ 
ever, the SBF relies on high aspect ratios to accurately 
model elongated particles, and other boundary conditions 
than periodic or free space are not as straight forward as 
in the LBM. The LBM can in contrast simulate complex 
geometries and nearly arbitrary particle shapes, only re¬ 
stricted by the pe in this respect. Moreover, the LBM is 
able to consider inertia, which allows physically more real¬ 
istic simulations and a better comparison to experiments. 
The high parallel efficiency of the LBM allows its execu¬ 
tion on massively parallel clusters and compensates for its 
higher computational effort. 

A possible extension of the SBF is the implementation of 
wall boundary conditions. Wall treatments in a boundary 
integral setting have been found to be challenging, and 
the literature is sparse. The most straightforward way to 
include these outer boundaries is to treat the wall in the 
same way as the immersed objects. A boundary integral 
over the wall is incorporated in the formulation and is 
discretized using special quadrature. This approach yields 
extra unknowns that must be computed and thus increases 
computational cost. The number of additional unknowns 
will depend on the size of the domain as well as on the 
resolution. A different way of including wall boundary 
conditions in the SBF is to use the method of images [SO]. 
However, this method is only feasible for one or two par¬ 
allel plane boundaries. 

In the future we may compare the two examined meth¬ 
ods for elongated particles in flows subject to wall bound¬ 
ary conditions and simulate the interactions of many sed¬ 
imenting particles, including LBM simulations at higher 
Reynolds number. The LBM may also be applied to sim¬ 
ulate collisions of elongated particles with walls and their 
motion in complex geometries. Further work could include 
multiphysics simulations of charged elongated particles in 
fluids subject to electric fields and their deposition on a 
charged surface, similar to the simulations with spheres 
in PO] . 


Appendix A. Detailed LBM results for single par¬ 
ticle motion 

Additional details on LBM results from Sec.|5] 

Table A.7. LBM sedimentation velocities of spherocylinders 

with different aspect ratios 1/e and radius r = 4dx (dx = 10“^ m), 
in [2560 dx]^ x 2688 dx domain with free-slip boundaries, and angular 
velocities ^lbm [816 dx]^ domain with no-slip boundaries. Fluc¬ 
tuations due to obstacle mappings (5u = ~ ^min)/^LBM 

6cj (defined analogeously) are shown, together with Reynolds num¬ 
bers based on f^LBM diameter (Rep^d) or length (Rep^^), or tip 
velocity (Rctip^d^ LBM velocities for lengthwise/sidewise 

translational- and rotational motion are compared to analytical so¬ 
lutions by Tirado et al. m for cylinders with lengths including (indi¬ 
cated by ‘Tir,wC’) and excluding (‘Tir,nC’) spherocylinder end-caps. 
Differences between LBM and Tirado solutions are shown as relative 
deviations A^f/ and A^uj w.r.t. the Tirado solutions. 
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